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The activity of natural killer (NK) cells is regulated by various processes including edu- 
cation/licensing, priming, integration of positive and negative signals through an array 
of activating and inhibitory receptors, and the development of memory-like functionality. 
These processes are often very complex due to the large number of different receptors 
and signaling pathways involved. Understanding these complex mechanisms is therefore 
a challenge, but is critical for understanding NK cell regulation. Mathematical approaches 
can facilitate the analysis and understanding of complex systems. Therefore, they may 
be instrumental for studies in NK cell biology. Here we provide a review of the different 
mathematical approaches to the analysis of NK cell signal integration, activation, prolifer- 
ation, and the acquisition of inhibitory receptors. These studies show how mathematical 
methods can aid the analysis of NK cell regulation. 
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INTRODUCTION 

Natural killer (NK) cells are innate lymphocytes that are impor- 
tant for early immune responses against viral pathogens and 
transformed cells (Vivier etal., 2008). They can be stimulated by 
cytokines or by direct contact with infected or transformed cells, 
resulting in cytolysis of the recognized ceO and the production of 
immune-modulatory cytokines. Target cell recognition is medi- 
ated by an array of stimulatory and inhibitory NK cell receptors 
(Lanier, 2008). The opposing signals induced by these receptors 
regulate NK cell activities (Watzl and Long, 2010). Inhibitory 
receptors are necessary to ensure NK cell self-tolerance in two 
ways. During NK cell development inhibitory receptors specific 
for self-MHC class I are necessary for the generation of functional 
competent NK cells in a process termed "education" (Elliott and 
Yokoyama, 2011). In mature NK cells the presence of inhibitory 
receptors prevents the attack of healthy autologous cells (Karre, 
2008). To study the acquisition of inhibitory receptors, their 
exact role during NK cell education and their interference with 
activating signaling is essential for the understanding of NK cell 
regulation. 

The immune system is one of the most complex biological sys- 
tems known. It involves various cell types including lymphocytes, 
each of which can exhibit a diverse repertoire of individual cells: 
T and B cells express clonally distributed somatically rearranged 
antigen receptors, whereas each NK cell can express different 
inhibitory and activating receptors. These cells do not act as 
individuals, but interact with other immune cells in a com- 
plex inter-ceUular network. Furthermore, lymphocyte activation 
is regulated by a complicated intra-cellular signaling machin- 
ery. As a result, the immune system and its functionality are 
highly complex and non-linear. The analysis of immune system 
data is thus non-trivial and calls for the use of systems biology 



methods, including mathematical and computational models 
(Walzer etal, 2007). 

MATHEMATICAL METHODS TO ANALYZE AND MODEL 
IMMUNE SYSTEM DATA 

Obtaining, analyzing, and comparing data on lymphocyte reper- 
toires yields many new insights (reviewed in Mehr etal., 2012). 
However, the formidable challenge is to formulate a description 
of the system in question on several levels: genetic, molecular, cel- 
lular, and systemic. Of course, no model can include all aspects 
of the real system - nor should it attempt to, because models 
must be simple enough to enable an exact analysis. It is up to the 
researcher to decide which aspects of the system to focus on, and 
which aspects to neglect (Mehr, 2001, 2005). 

Before performing any multi-scale modeling, one should know 
how to model events on a single scale. These methods can be 
combined in multi-scale modeling in various ways, using for 
each scale the most appropriate method (Shannon and Mehr, 
1999; Kirschner etal, 2007; Materi and Wisliart, 2007). The 
main modeling techniques used so far in immunology are: mod- 
els of population dynamics using ordinary differential equations 
(ODEs), including delay differential equations (DDE), or - in 
case the numbers of objects (e.g., cells or molecules) is not very 
large and stochastic effects are important - stochastic differential 
equations (SDEs). The latter (and other so-called "Monte Carlo" 
modeling methods, such as stochastic agent-based models - see 
below) have been developed to better account for intrinsic cel- 
lular and molecular fluctuations. For example, see our (mostly 
ODE) models of lymphocyte development - T cell development 
in the thymus (Mehr etal, 1993, 1994, 1995, 1996a,b, 1997, 1998) 
and homeostasis in peripheral blood (Mehr and Perelson, 1997), 
and ODE or hybrid models of B cell development and responses 
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(Mehr et al, 2003; Shahaf et al., 2004, 2006, 2010). Stochastic mod- 
eling has also been extensively used by several groups, e.g., the 
Hodgkin group's "Cyton" model (Lee et al., 2009; Markham et al, 
2010; Duffy and Hodgkin, 2012; Duffy etal, 2012). In case spa- 
tial effects are important, one can either divide the space into 
compartments and stiU use ODE modeling for each compart- 
ment, or use partial differential equations (PDE), which explicitly 
model the dependence of the model's entity on space. A differ- 
ent method which has become popular with the availability of 
high-speed computation, and is also more intuitive and accessible 
for biologists, is agent-based modeling (ABM) - where individ- 
ual entities (agents) such as molecules, cells or individuals are 
explicitly modeled as program objects, which move, respond to 
their environment and/or change their internal state based on 
pre-defined rules (Chavali etal., 2008). Stochasticity is easy to 
buUd into ABM, as are different levels of description. One can 
also create hybrid models where different levels of description 
are modeled using different model types (Shahaf etal., 2008). 
Many such models of lymphocyte migration or molecular pro- 
cessed were published (e.g., Beltman etal., 2011; Kaplan etal, 
2011). A review of ODE, DDE, PDE, stochastic, and ABM mod- 
els in immunology, including a comparison of their relative 
advantages and disadvantages, has recently been published (Kim 
etal., 2009). 

MODELING NK CELL SIGNALING PATHWAYS AND THE 
DYNAMICS OF NK CELL ACTIVATION 

The understanding of T cell receptor signaling has been greatly 
aided by mathematical approaches (for examples, see Feinerman 
et al., 2008; Das et al, 2009). NK cells rely on similar signaling path- 
ways for their activation. However, NK cell activation is mediated 
by several different activating receptors, which can synergize with 
each other (Bryceson etal, 2006). Additionally, inhibitory recep- 
tors can counteract NK cell activation (Watzl and Long, 2010). 
This makes NK cell regulation even more complex. Analyzing and 
predicting the behavior of NK cell repertoires necessitates models 
that take into account the complexity of whole NK cell reper- 
toires (Johansson et al, 2009), as well as the functional differences 
between NK cell repertoire subsets (Brodin et al., 2010, 2012) and 
the molecular processes underlying these differences. Mathemati- 
cal modeling may therefore help to understand the integration of 
these opposing signals. However, the lack of detailed knowledge 
about the complex topology of NK cell signaling pathways presents 
a challenge for modeling approaches. 

In a purely theoretical approach using stochastic agent-based 
simulation one study has developed a detailed molecular model of 
NK cell activation, incorporating membrane-proximal signals and 
affinities of receptor-ligand interactions (Das, 2010). The engage- 
ment of activating receptors results in the activation of the Vav 
and Erk signaling pathways. Interestingly, the study showed that 
on a single cell level the responses can be digital in nature so that 
the latter signaling pathways are either active or not. Such digi- 
tal behavior of Erk phosphorylation matches previous results in 
T cells (Altan-Bonnet and Germain, 2005), where opposing feed- 
back loops of positive (Erk) and negative (SHP- 1 ) signals created a 
digital response. However, the model showed that in NK cells such 
a digital response is possible even in the absence of any feedback 



loop (Das, 2010), because of the presence of inhibitory receptors 
in the NK cell synapse. 

Another study used an ODE-based "ensemble modeling" 
approach (Kuepfer et al., 2007) combined with experimental veri- 
fication to investigate the interplay of activating and inhibitory sig- 
nals in NK cells (Mesecke etal., 2011). This study experimentally 
confirmed a digital Vav phosphorylation and de-phosphorylation 
induced by activating and inhibitory signaling, respectively. The 
latter result suggests that Vav phosphorylation is an integration 
point at which activating and inhibitory signaling converge (Steb- 
bins et al., 2003). The digital nature of Vav phosphorylation (Das, 
2010; Mesecke et al, 201 1) could explain how these opposing sig- 
nals are integrated to yield a yes or no decision about NK cell 
activity, resulting in the survival or the death of a conjugated 
target cell. 

The inhibition of NK cells depends on the amount of MHC 
class I expressed by the target cell. It has been shown that there is 
a clear threshold in the amount of surface MHC class I required 
for inhibition of NK cell cytotoxicity. A study by Almeida et al. 
modeled the dependence of the activation/inhibition threshold 
on MHC class I and inhibitory receptor expression levels at the 
single cell level using ODEs (Almeida et al, 20 1 1 ) . They fitted their 
model to data on NK cell killing of target cells with various MHC 
levels, obtained in experiments from many different human NK 
cell clones, in order to see whether clones differ in their activation 
threshold, maximal killing capacity, or both. The results showed 
that NK cells differ mostly in their intrinsic activation threshold, 
which may be set during the process of NK cell education. 

An agent-based modeling study by Kaplan etal. (2011) 
addressed the question how signals from activating and inhibitory 
receptors are spatially integrated within the NK cell immunologi- 
cal synapse. The study constructed a computer simulation of the 
assembly of the NK cell-target cell synapse, which was based on 
many experimental findings. The NK cell immune synapse was 
modeled as two parallel two-dimensional grids, representing the 
contact areas on the membranes of the two interacting cells - the 
NK and the target cell grids. The simulation suggested that the 
NK cell does not simply sum up all activating and inhibitory sig- 
nals over the entire synapse. Instead, a model considering that 
the inhibitory signal is locally confined to a radius of 3-10 recep- 
tor molecules was closer to the experimentally observed behavior. 
This is in line with the fact that inhibitory receptors mediate their 
negative signal through direct recruitment of a protein tyrosine 
phosphatase such as SHP-1, and suggests that the active phos- 
phatase cannot diffuse very fast. Therefore, inhibitory receptors 
have a locally limited inhibition radius. 

The simulation by Kaplan et al. (201 1) further showed that the 
concerted motion of receptors in clusters - such as those formed by 
molecule aggregation on certain membrane microdomains - sig- 
nificantly accelerates the maturation of the NK cell immunological 
synapse. Membrane microdomains have also been implicated by 
the study of Mesecke et al. The mathematical modeling showed 
that the association of Src-famUy kinases with activating but not 
with inhibitory receptors was essential for the experimentally 
observed Vav activation (Mesecke etal., 2011). While this might 
suggest a direct association between activating NK cell receptors 
and Src-family kinases, it could also support the known role of 
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membrane microdomains in NK cell activation (Lou et al., 2000; 
Fassett etal, 2001; Watzl and Long, 2003; Endt etal, 2007; Guia 
etal, 2011). As these domains are enriched in Src-family kinases, 
the recruitment of activating receptors to these domains might 
be the event that brings kinases and receptors in close proximity. 
Additionally, clustering of activating receptors would make them 
more resistant to the locally restricted negative signal of inhibitory 
receptors (Kaplan etal., 2011). 

A modeling study by Burroughs etal. (2011) addressed the 
question of receptor segregation at the NK cell immunological 
synapse. They specifically asked the question of whether the seg- 
regation of different receptors at the immunological synapse is 
simply a result of the different sizes of the proteins. They used 
a novel energy model that can analyze protein redistribution 
and developed a methodology for analyzing and quantifying data 
from two color fluorescence images. Using these methods, they 
showed that the observed protein segregation patterns in actual 
synapses can be explained by differences in protein size alone. 
This was independent of the actin cytoskeleton and could explain 
how certain receptor patterns form within the immunological 
synapse. 

Once a NK cell has been sufficiently stimulated by a locally 
attached target cell, a cytotoxic process is initiated, resulting in the 
lysis of the target cell. A study by Almeida et al. (20 11 ) , explored the 
kinetics of NK cell activation and target cell lysis at the population 
level. Experiments have shown that while the number of NK cell- 
target cell conjugates reaches saturation within 30 min, significant 
target cell lysis is only observed after 3-5 h. A model assuming 
that the NK cell can kill its target cell only after being primed by a 
former target cell encounter could reproduce the experimentally 
observed kinetics and explain this delay. Thus, the population of 
NK cells is divided into resting NK cells that do not kill targets, 
and activated NK cells that are capable of killing. This suggests 
that NK cells need to be primed before they can kill target cells, 
and that this priming occurs after the first encounter of the resting 
NK cell with a target cell. 

NK CELL POPULATION DYNAMICS 

Mathematical modeUng has been used in several studies to eluci- 
date the dynamics of NK cell populations, their development and 
activation. Early mechanistic modeling studies of NK cell pop- 
ulation dynamics (as opposed to simple curve-fitting studies of 
killing dynamics) included attempts to clarify their role in immune 
surveillance (Merrill, 1981, 1983) and in immunotherapy (Nani 
and Oguztoreli, 1994), or elucidate the dynamics of their adhesion 
to target cells (Kuznetsov, 1996). In recent years, there have been 
several more mathematical modeling studies on the role of NK 
cells in response to infections (Hancioglu et al., 2007; Wodarz et al, 
2007; Elemans et al, 2011) and cancers (Cappuccio et al., 2006). In 
many of these studies, however, the models focus on the adaptive 
T cell response, and NK cells are just an additional feature. We 
are only aware of one study, which focused primarily on NK cell 
population dynamics and turnover (Zhang etal., 2007). Recently, 
one report investigated the proliferation of NK cells in response to 
IL-15 (Zhao and French, 2012). IL-15 plays an important role in 
NK cell development, survival and homeostasis, but also during 
virus-induced NK cell proliferation. The authors developed a 



series of two-compartment models describing the dividing and 
non-dividing NK cell subpopulations. NK cells could be recruited 
from one compartment into the other and possessed distinct pro- 
liferation and cell death rates within a compartment. Comparing 
experimental results to different model variations showed that IL- 
15 modulates the recruitment rate from the non-dividing to the 
dividing compartment and influences the division rates. However, 
it did not have any influence on the time until the first division. 
The model could even predict NK cell proliferation dynamics in 
an environment of changing IL-15 concentrations, as it would be 
the case during viral infections, demonstrating the usefulness of 
such mathematical models. 

NK CELL REPERTOIRE FORMATION 

The dynamics of NK cell repertoire development and education 
are even more complex than the population dynamics of NK cell 
subsets. Initially, two conceptual models have been proposed by 
Vance and Raulet (1998), to explain the process of NK cell "educa- 
tion" in which the cells adapt to the self-MHC environment. Both 
models have assumed that in order to mature, there should be at 
least one self-MHC-specific receptor on the NK cell (the "at least 
one" hypothesis). The first is the sequential activation model. In 
this scheme, genes for Ly49 inhibitory NK cell receptors are ran- 
domly expressed; once a receptor gene is expressed, it is known 
to remain stably activated. During maturation, the cell expresses 
new receptors, until a receptor that binds self-MHC is expressed. 
The cells are periodically tested for interaction with self-MHC 
molecules. Strong interactions between the receptors and their 
MHC ligands prevent additional receptors from being expressed 
and result in the maturation of the cell. In this scheme, a sin- 
gle testing step accomplishes the education process; however this 
single step may be repeated many times during the cell's devel- 
opment. The second scheme is the two-step selection model, in 
which the repertoire is fully formed at an initial stage by a stochas- 
tic process of receptor expression, and subsequently shaped by 
two selection steps: one selects for cells expressing at least one 
self-specific receptor, and the other selects against cells expressing 
multiple self-specific receptors. The two-step selection thus occurs 
only once for each cell, when it has completed its receptor gene 
activation. Depending on the signaling thresholds at these steps, 
the process may allow maturation of cells expressing more than 
one self-specific receptor. 

Mathematical methods are ideally suited to distinguish between 
these two models. In the study by Salmon-Divon etal. (2003a), the 
"two-step" model was implemented as a computer program which 
calculates the expected frequencies of cells with given receptor 
expression patterns as function of the selection thresholds. The 
"sequential" model of NK cell selection was implemented as a 
stochastic agent-based simulation of the development of the NK 
cell repertoire, again giving the expected frequencies of cells with 
given receptor expression patterns based on model parameters. 
The original Raulet group data was not sufficient to conclu- 
sively decide between the models (Salmon-Divon etal., 2003b). 
However, the use of specifically generated, larger data set from 
single-MHC mice showed that the two-step selection model fits 
the data significantly better than the sequential model (Johansson 
et aL,2009). 



www.frontiersin.org 



December 2012 | Volume 3 | Article 359 | 3 



Watzl etal. 



Mathematical approaches to NK cell regulation 



Additional studies aimed to determine whether the MHC class 
I background affects NK cell repertoire composition, in both 
humans and mice. For this purpose, data on NK cell repertoires 
extracted from the bone marrow (BM) of mice from different 
single-MHC and MHC-deficient backgrounds was analyzed. In 
addition, NK cell repertoire data from blood samples of humans 
with different MHC backgrounds was used. In both cases, cells 
were probed for the expression of five inhibitory receptors (from 
the Ly49 or KIR families, in mice and humans, respectively). The 
results of these studies are described in another paper in this 
issue (Simon, M. etal, submitted). 

OUTLOOK 

There are still many challenging and open questions about differ- 
ent aspects of NK cell regulation to date. We still need to learn more 
about NK cell education, subsets, signaling and even memory 
(Paust and von Andrian, 2011). The studies described above have 
demonstrated that mathematical approaches can help to answer 
some of these open questions. They are ideally suited to analyze 
complex data sets, make predictions about complicated systems 
or simply aid experimental design. We therefore expect to see 



many more studies involving mathematical approaches to NK cell 
regulation in the near future. 
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